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Abstract 

We consider the problem of pricing path-dependent options on a basket of underlying 
assets using simulations. As an example we develop our studies using Asian options. 

Asian options are derivative contracts in which the underlying variable is the average 
price of given assets sampled over a period of time. Due to this structure, Asian op- 
tions display a lower volatility and are therefore cheaper than their standard European 
counterparts. 

This paper is a survey of some recent enhancements to improve efficiency when 
pricing Asian options by Monte Carlo simulation in the Black-Scholes model. We 
analyze the dynamics with constant and time-dependent volatilities of the underlying 
asset returns. 

We present a comparison between the precision of the standard Monte Carlo method 
(MC) and the stratified Latin Hypercube Sampling (LHS). In particular, we discuss the 
use of low-discrepancy sequences, also known as Quasi- Monte Carlo method (QMC), 
and a randomized version of these sequences, known as Randomized Quasi Monte Carlo 
(RQMC). The latter has proven to be a useful variance reduction technique for both 
problems of up to 20 dimensions and for very high dimensions. 

Moreover, we present and test a new path generation approach based on a Kronecker 
product approximation (KPA) in the case of time-dependent volatilities. KPA proves 
to be a fast generation technique and reduces the computational cost of the simulation 
procedure. 

Key Words: Monte Carlo and Quasi-Monte Carlo simulations. Effective dimensions. 
Path-generation techniques. Path-dependent options. 

1 Introduction 



The financial industry has developed a variety of derivative contracts in order to fulfil 
different investor needs. Path-dependent options play a fundamental role in financial 
engineering and can display different exotic features. 

Exotic contracts that are widely used are Asian options, barrier options and look- 
back options both with American and European style. An un-biased and efficient 
pricing procedure is fundamental and a vast research is done in order to obtain fast 
and efficient estimations. Common approaches rely on finite differences methods and 
Monte Carlo simulations. 

Finite differences methods consist in discretizing the partial differential equation 
whose solution gives the price of the options while Monte Carlo methods face the 



problem from a probabilistic point of view. It estimates the price as an expected value 
by its integral formulation. 

The former method returns the fair price of the option for different times and val- 
ues of the underlying variable but is practically unfeasible for complicated multi-asset 
dependence. 

On the other hand, Monte Carlo simulation calculates the fair price in a single time 
point and can be applied to various situations. 

Its fundamental property is that its order of convergence is 0{1/ ^/n) and does not 
depend on the number of random sources of the problem. Although it does not display a 
high order of convergence, it proves to be efficient for pricing complex exotic contracts. 

The aim of this report is to describe standard and advanced Monte Carlo techniques 
applied for multi-asset Asian options of European style. In particular we concentrate 
our studies to stratification and Quasi-Monte Carlo approaches. 

Standard Monte Carlo can be seen as a numerical procedure aimed to estimate 
integrals in the hypercube [0, 1]*^ by generating different scenarios with uniform random 
variables. Stratification achieves the same task by drawing uniform random variates in 
a smaller set in [0,1]'' introducing correlation. 

Quasi-Monte Carlo methods drop off all probabilistic considerations and focus on 
the problem of generating a sequence of points that uniformly covers the hypercube 
[0,1)'* (the theory is built-up for right-opened intervals). The sequence is absolutely 
deterministic and different drawings lead to the same points. 

From the mathematical point of view, it introduces the concept of discrepancy and 
star-discrepancy that quantify how well the sequences cover [0, 1)''. 

Hawkla and Koksma proved the fundamental inequality, named after them, that 
provides the bound for the estimation error of the target integral depending on the 
discrepancy. 

Low-discrepancy sequences are those whose estimation error is O(^). The con- 
vergence rate depends on the dimension d and is lower than the Monte Carlo error for 
small d. There exist several low discrepancy sequences, among them are the Halton, 
the Faure, the Sobol and the Niederreiter-Xing sequences. A fundamental reference on 
this topic is Niederreiter [T7| . 

Quasi-Monte Carlo methods can be unpractical because the computation of the 
error is potentially more difficult than the estimation of the target integral while the 
better uniformity can be lost even for low values of d. 

Standard Monte Carlo, stratification and Quasi-Monte Carlo methods form a hier- 
archy for the generation of uniform points. 

A further step ahead can be taken by randomizing these sequences while preserv- 
ing the low-discrepancy. This technique is called scrambling, Owen [21] provides an 
extensive description on the subject. 

The application to options pricing is straightforward. Standard models for price 
dynamics involve multidimensional Ito processes so that pricing exotic contracts might 
require a high-dimensional integration. It necessitates careful implementation of the 
simulation especially when Quasi-Monte Carlo methods are used. 

Many works have been done to investigate the problem. Acworth, Broadie, and 
Glasserman [T] provided a first comparisons between variance reduction techniques and 
Quasi-Monte Carlo methods and Caflisch, Morokoff and Owen [3] analyzed the effec- 
tive dimension of the integration problem for mortgage-backed securities by ANOVA 
considerations. Caflisch, Morokoff and Owen [i] and Owen [12] showed that only few 
random sources really matter and suggested to choose for them a better generation 
technique. 

We focus our investigation on pricing Asian options in a multi-dimensional Black- 
Scholes model both for constant and time-dependent volatilities. In this framework, 
standard path-generating techniques are the Cholesky decomposition, the principal 
component analysis (PCA) and linear transform LT. The last two have been proven 
to be essential for ANOVA in order to recognize effective dimensions so that an efficient 
RQMC can be run. 
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When constant volatilities are considered, the path-generation procedure can be 
simplified relying on the properties on the Kronecker product while this is not possible 
for time-depending volatilities. 

As for this task, we propose a new approach based on a Kronecker product approxi- 
mation. The general problem consists in approximating the global correlation matrix of 
the price returns into the Kronecker product of two smaller matrices. We assume that 
the former of the two is the auto-covariance matrix of a single brownian motion. Indeed, 
we suppose that most of the variance of the global process is carried out by each driving 
brownian motion. The latter matrix would be an approximation of the total covariance 
matrix among the asset returns during the lifetime of the contract. The original and 
target path is reobtained by Cholesky decomposition. As for this last step we develop 
an ad hoc realization of the Cholesky decomposition suited for the global correlation 
matrix. This procedure is intended to reduce the computational burden required to 
evaluate the whole set of eigenvalues and eigenvectors of the global covariance matrix. 

The last step of the simulation is the computation of the Asian price via simulation 
using standard Monte Carlo, LHS and RQMC approaches. As for the last one we 
perform a Faure-Tezuka scrabling version of the Sobol' sequence, which is the most 
used low discrepancy sequence in finance. 

In the case of constant volatility, we set our investigation as in Dahl and Benth , 
[7] and Imai and Tan [TT]. We compare our results and analyze the precision of the 
simulation for different path-generation methods and Monte Carlo approach. 

As for the time-dependent volatility market we test the KPA method and compare 
its results with those obtained with the PCA decomposition. 

Summarizing the report evolves as follows: section 2 introduces the market. Section 
3 describes the pay-off of Asian options and presents the problem as an integral for- 
mulation. Section 4 defines effective dimensions in truncation and superposition sense. 
Section 5 defines the Kronecker product and list some of the main properties. Section 
6 describes some path-generation procedures and in particular, it introduces the KPA 
method. Section 7 is a brief introduction to low-discrepancy sequences and scrambling 
techniques. Section 8 describes the simulation procedure we adopt. Section 9 shows 
and comments the estimated results for different scenarios both in the constant and 
time-dependent cases. 

2 The Market 

We consider a complete, standard financial market 9Jt in a Black-Scholes framework, 
with constant risk-free rate r and time-dependent volatilities. There are M + 1 assets 
in the market, one risk free asset and M risky assets. The price processes of the assets 
in this market are driven by a set of stochastic differential equations. 

Suppose we have already applied the Girsanov theorem and found the (unique) risk- 
neutral probability, the model for the risky assets is the so called multi-dimensional 
geometric brownian motion: 

So[t) = e'-* (1) 

dS,{t)^rS,{t)dt + a,{t)S,{t) dW,{t), i^l,...,M. (2) 

Here Si (t) denotes the i-th asset price at time i, ai (t) represents the instantaneous time- 
dependent volatility of the z-th asset return, r is the continuously compounded risk-free 
interest rate, and W (t) — {Wi (i) , . . . , Wm {t)) is an M-dimensional Brownian motion. 
Time t can vary in R!j_, that is, we can consider any maturity T G M.'^ for all financial 
contracts. 

The multi-dimensional brownian motion W (t) is a martingale, that is, each com- 
ponent is a martingale, and satisfies the following properties: 

E[W,{t)]=0, z = l,...,M. 

[W^,Wk]{t) = p,kt, z,fc = l,...,M. 
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where represents the quadratic variation up to time t and pik the constant 

instantaneous correlation between Wi and Wk ■ 

Consider a generic maturity T, we can define a time grid T = {<i, . . . ,iAr} of N 
points such that ti < t2 < ■ ■ ■ ,tN = T, we recall that the sampled covariance matrix 
Ri^m = E [Wi (ti) Wi (tm)], l,m — 1, . . . , N oi each Brownian motion in equation ^ is: 

/ h h ... h \ 

ti t2 ... t2 



R 



(3) 



\ ti t2 ... In / 

This matrix is symmetric and its elements Ri^m = U /^tm have the peculiarity to be 
constant after reflection about the diagonal. We will refer to this feature as boomerang 
shape property. 

In order to complete the picture of our environment, we need to define the matrix 
whose elements are T^i^kit) =Pik<^iit)o'k(t), i,k — 1,... M. This is a time de- 
pendent covariance matrix evolving according to the dynamics of the time-dependent 
volatilities and the constant correlation among the asset returns. 

Avoiding all the calculation (see Rebonato [21] and Glassermann ^ for further 
details), we derive the global covariance matrix Em at that assumes the expression below: 



E(t2 



(4) 



The global covariance matrix is very simple and enjoys the boomerang shape property 
with respect to the block-matrix notation. All the information is carried out by N 
time-varying M x M matrices. 

Each element depends on four indexes: 

[J:mn),i.]_ = / ai{t)ak{t)pikdt (5) 



with i,k = 1, . . . , M and l,m = 1, . . . ,N. 

Applying the risk-neutral pricing formula, the value at time t of any European 
T-maturing derivative contract is: 

V{t)^exp{r{T-t))E[cj){T)\Tt]- (6) 

E denotes the expectation under the risk neutral probability measure and 4'{T) is a 
generic J-t measurable function, with Tt = (j{0 < t < T;W{t)}, that determines the 
payoff of the contract. Although not explicitly written, the function (t>(T) depends on 
the entire multi-dimensional brownian path up to time T. 



3 Problem Settlement 

We will restrict our analysis to Asian options that are exotic derivative contracts that 
can be written both on a single security and on a basket of underlying securities. 
Hereafter we will consider European-style Asian options whose underlying securities 
coincide with the M + 1 assets on the market. This is the most general case we can 
tackle in the market 9Jl, because it is complete in the sense that we can hedge any 
financial instrument by finding a portfolio that is a combination of this M + 1 assets. 



3.1 Asian Options Payoff 

The theoretical definition of Asian options price is: 



a^{t) ^ exp{r{T -t))E 



lo S^{t)dt 



T 



-K 



Option on a Single Asset (7) 
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a{t) = exp{r{T -t))E 



lo Y.f=iwA{t)dt 



T 



K 



Option on a Basket 



(8) 



K represents the strike 
are usually 



where we assume that the start date of the contract is t = 0. 
price and coefficients Wi satisfy X^fii '^i — 1- Contingent claims ^ and 
referred as weighted Asian options. 

In practice no contract is agreed according to equations ([7]) and The integrals 
are approximated by sums; often these approximations are written in the contracts by 
specifying the number and the sampling points of the path. 

Approximation for ([7|) and ([8]) can be carried out by using the following expressions: 



Qi (t) = exp {r{T - t)) E 



a{t) = exp (r{T - t))'. 



where coefficients Wij satisfy 



N 



K 




= 1. 



Tt 



Option on a Single Asset 

(9) 

Option on a Basket (10) 



European options with payoff functions (O and (fTU)) are called arithmetic weighted 
average options or simply arithmetic Asian options. When M > and TV = 1 the 
payoff only depends on the terminal price of the basket of M underlying assets and the 
option is known as basket option. 

No closed-form solution exists for Asian options arbitrage-free price, neither for 
single nor for basket options both for theoretical and finitely monitored payoff. In 
order to obtain a correct valuation of the price we are compelled to turn to numerical 
procedures such as the Monte Carlo estimation or the finite difference methods. 

The latter is based on a convenient and correct discretization of the partial differ- 
ential equation associated to the risk neutral pricing formula via the Feynmann-Kac 
representation. The finite difference method returns the price for all the times and 
initial values of the underlying assets. Vecer [55] and [27] found a convenient approach 
for the single asset case and presents the comparison with other techniques. The main 
drawback is the stability of the method that is practically unfeasible for options on a 
basket. 

Monte Carlo simulation is a numerically intensive methodology that provides unbi- 
ased estimates with convergence rate not depending on the dimension of the problem 
(the number of random sources to draw) . The cases of high values for the problem di- 
mension find interesting applications in finance including the pricing of high-dimensional 
multi-factor path-dependent options. In contrast to the finite difference technique, the 
Monte Carlo method returns the estimate for a single point in time. It is a fiexible ap- 
proach but requires ad hoc implementation and refinements, such as variance reduction 
techniques, in order to improve its efficiency. 

The main purpose of the standard Monte Carlo method is to numerically estimate 
the integral below: 



/ = 



[0,1]' 



/(x) dx. 



(11) 



The integral / can be regarded as E [/(?/)], the expected value of a function /(.) of the 
random vector U that is uniformly distributed in hypercube [0, 1]''. 

Monte Carlo methods simply estimate / by drawing a sample of n independent 
replicates Ui . . . ,Un of U and then computing the arithmetic average: 



^ ^ 1 " 



(12) 



The Law of Large Numbers ensures that /„ converges to / in probability a.s. and 
the Central Limit Theorem states that / — /„ converges in distribution to a normal with 
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mean and standard deviation a / \/n witli = y Jq (/(x) — dx. The convergence 
rate is than 0{l/^/n) for all dimensions d. The parameter a is generally unknown in 
a setting in which / is unknown, but it can be estimated using the sampled standard 
deviation or root mean square error (RMSE): 



1 2 

— E(/(f^»)-^«) ■ (13) 

i=l 

Refinements in Monte Carlo methods consist in finding techniques whose aim is to 
reduce the RMSE, known as variance reduction techniques, without changing the con- 
vergence rate. In contrast, the Quasi Monte Carlo version focuses on the improvement 
of the convergence rate by generating sequences in [0, 1]'' with high stratification in 
order to uniformly cover the hypercube. These sequences are no longer random and 
estimates and errors are not based on probabilistic considerations. 

As far as our case is concerned, we need to formulate the problems (O and (fT(]| for 
pricing Asian options as integrals of the form pip in order to apply the Monte Carlo 
procedure. 



RMSE : 



\ 



3.2 Problem Formulation as an Integral 



The model 371, presented in the first section, consists of the risk-free money market 
account and M assets driven M geometric brownian motion described by equation ^ 
whose solution is: 



Si (t) = Si (0) exp 



ds 



a, (s) dW, (s) 



,M. 



(14) 



The quantity jj" ^^^ds is the total volatility for the i-th asset. The solution is 



a multi-dimensional geometric brownian motion, written GBMyr, "^'^ ds j , in the 
sense that it can be obtained applying Ito's lemma to Si{t) = f {Xi{t)) = 



ith 



Xi[t) the i-th component of the multi-dimensional brownian motion with drift r and 



■j-th diffusion ^^^ds, written BM^r, J* ^^-^ds 



Under the assumption of constant volatility the solution is still a multi-dimensional 
geometric brownian motion with the following form: 



Si (t) = St (0) exp 



2 



t + <T,W^ (t) 



1, 



(15) 



In compacted notation the solution (|15p is GBM^^r, ^ 

Pricing Asian option requires to monitor the solutions and at a finite set 
of points in time {ti, . . . , ijv}- This sampling procedure yields the following expressions 
for time-dependent and constant volatilities: 



Si{tj) = Si{0)exp 



(16) 



Si{t-j) = St{0)exp 



(17) 



where the components of the vector (Zi(ti), . . . Zi{tN), ^2(^1), . . . , ZM{tN)) are M x N 
normal random variables with zero mean vector and covariance matrix Smw, whose 
form simplifies in the case of constant volatilities as we will be shown in Section 4. 
The payoff at maturity T of the arithmetic average Asian option is then: 



Pa{T) - (.g(Z) - K) 



(18) 
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where 



MxN 



g{Z) = ^ exp {p.k + Zk) 



(19) 



and 



^fc = \T^{wklk2Sk^{Q)) + 



( 



r 




(20) 



for constant volatilities or 



/ife = ln(wfejfe2S'fci (0)) + rtk^ 



(21) 
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for time-dependent volatilities. The indexes ki and fc2 are ki = (k — l)modM, k2 = 
[(fc— 1)/M]-|-1, respectively, where mod denotes the modulus and [.] the greatest integer 
less than or equal to x. 

The calculation of the price a[t) in equation ()10l) can be formulated as an integral 
on 'U.^^ in the following way (see Dahl and Benth [6j and f^): 



Fz is the cumulative distribution of the normal random vector -/V(0, S^vm)- 

In the following section we will show how to obtain the random vector Z starting 
from a vector of independent and normally distributed random variables e. Once this 
generation is carried out, we can apply the inverse transform method to formulate the 
pricing problem as an integral of uniform random variables in the hypercube [0, l]^^^ 
and use Monte Carlo methods: 



In the following sections we will present recent enhancement based on ANOVA for 
high dimensional Monte Carlo and Quasi-Monte Carlo simulations in order to estimate 
the integral for the pricing Asian option on a basket of underlying assets both for 
constant and time-dependent volatilities. 

4 Effective Dimensions 

When the nominal dimension d of the problem of estimating the integral (jlip is one, 
there are standard numerical techniques that give a good accuracy when / is smooth. 
Considerable problems arise when d is high. 

Recent studies proved that many financial experiments present problem dimensions 
lower than the nominal one. Owen (1998) J_9J and Caflisch, Morokoff and Owen [3] 
studied the application of ANOVA for high-dimensional problems and introduced the 
definition of effective dimension. It is possible to study some mathematical properties 
of the function / and try to split it in order to reduce the computational effort. The 
ANOVA decomposition consists of finding a representation of / into orthogonal func- 
tions each of them depending only on a subset of the original variables. This is the 
peculiar and stronger condition that makes ANOVA different and more powerful with 
respect to the usual Least Squared method. 

Let A — {!,..., d} denote the set of the independent variables for / on [0,1]''. 
/ could be written into the sum of 2'^ orthogonal functions each of them defined in 
different subsets of A, that is depending only on the variables in each of these subsets: 




(22) 




(23) 




(24) 



u<ZA 
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Now let |u| denote the cardinality of u, Xu the |u|-tuple consisting of components Xj 
with j € u, and —u being the complement of u in A. Then set the function as: 

/«(x) = / (f{z) - fvi^)) dz u (25) 

Equation (|25p defines /„ by subtracting what can be attributed to the subsets of u, 
and then averaging over all components not in u. In the function setting /u(xu) only 
depends on Xu- 

Denoting cr^ = /(/(x) — /)^ dx, crj = //i,(x)-^dx, (Tq = 0, supposing tr < +00 and 
|w| > it follows: 

- E (26) 

Equation (j26p partitions the total variance into parts corresponding to each subset u C 
A. The /„ exhibits some nice properties: if j G u the line integral Jj^ ^-^ fu{'^) dxj — 

for any Xk with k ^ j, and if m 7^ w / fu{x)fv{x) dx = 0. 

Exploiting the ANOVA decomposition the definition of effective dimension can be 
given in the following ways: 

Definition 1. The effective dimension of f, in the superposition sense, is the smallest 
integer ds such that J2o<\u\<ds '^u ^ P'^^ ■ 

The value ds depends on the order in which the input variables are indexed. 
Definition 2. The effective dimension of f , in the truncation sense, is the smallest 
integer dr such that J2uC{i dt} '^u — P^'^ ■ 

< p < 1 is an arbitrary level; the usual choice is p = 99%. 

The definition of effective dimension in truncation sense reflects that for some in- 
tegrands only a small number of the inputs might really matter. The definition of 
effective dimension in superposition sense takes into account that for some integrands 
the inputs might influence the outcome through their joint action within small groups. 
Direct computation leads: ds < dx < d. 



5 The Kronecker Product 



The Black-Scholes model was originally built up under the hypothesis of constant 
volatilities for all the assets. If this assumption drops off the main ideas underlying 
the market Tl described above do not change and fundamental results still hold. The 
constant volatility case reduces the computational complexity of the analysis and sim- 
plifies many calculations. 

In the following we present some useful properties of the brownian motion, of its 
sampled autocovariance matrix and of the global covariance matrix. Furthermore, we 
introduce the Kronecker product that will prove to be a powerful tool for reducing the 
computation burden and a fast way to generate multi-dimensional brownian paths. 

The sampled covariance matrix of each brownian motion, R, enjoys many properties 
due to its particular boomerang form. We list some of them below: 

1. The inverse of i? is a symmetric tri-diagonal matrix: 



R- 



( 



tl(t2-tl} 



-tl 
-tl 



-*l)(t3 
1 



-t2) 



(t3- 



t3 
t4- 



-*2)(t4- 
1 

^4 — ^3 



■t3) 



^4 — ^3 



(t„-i 



-t„-2)(t„ 
1 

t„-t„_ 



(27) 
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is a sparse matrix and low memory is required to store it. R and share 
the same set of eigenvectors and have inverse eigenvalues (the matrices are both 
definite positive). 

2. The Cholesky decomposition of R gives a boomerang shaped matrix C. 

Definition 3 (Cholesky Decomposition). Given any hermitian, definite pos- 
itive matrix A, then A can be decomposed as: 



A = CaC*. 



(28) 



where Ca is a lower triangular matrix with strictly positive diagonal entries, and 
C* denotes the conjugate transpose of C. The Cholesky decomposition is unique 
and the Cholesky matrix can be interpreted as a sort of square root of R; as far 
as the Cholesky decomposition of a symmetric matrix A is concerned C% must 
be replaced by Cj. 

After direct computation Cr shows the form below: 







\ 



(29) 



y/tw — ijv-1 / 



In the case of an equally spaced time grid, the Cholesky matrix is just a lower 
triangular matrix whose elements are all equal to the time step At. 

3. The inverse of the Cholesky matrix is a sparse matrix, in particular it is a bi- 
diagonal matrix whose elements on the same row are equal and in opposite sign: 







V 















(30) 



All these results prove to be useful for the simulation and reduce the number of opera- 
tions for the brownian path generation. 

As for constant volatilities, both the covariance matrix among the asset returns and 
the global covariance matrix simplify and are not time-depending anymore. 

Let S be a covariance matrix depending on the correlation among the asset returns 
whose elements are: Ej^fc =PikO'icrk, i,k = 1, . . . M, then the global covariance matrix 
Smjv displays the following form: 



-iMN 



tiY, 

t2T. 



V tiE t2T. 



t2T, 
tN^ J 



(31) 



This matrix is obtained by repeating the constant block of covariance E at all the points 
of the time grid. 

This kind of mathematical operation is known as Kronecker product, denoted as (g). 
As such, Emjv can be identified as the Kronecker product between R and E, R^'E. 
The Kronecker product reduces the computational complexity by enabling operations 
on a {N X M,N X M) matrix using two smaller matrices that are N x N, and M x M 
respectively. 
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Definition 4 (The Kronecker Product). The Kronecker product of Am^xnA G 
^niAXTiA g^g^ BynsxiiB £ R™B><"B^ written A® B, is the tensor algebraic operation 
defined as: 



A<SiB 



I aiiB ai2B ... ai^^B ^ 
a2iB a22-B ... a2„^i? 

\ 0.mAlB dniAlB ... CLmAUAB J 



(32) 



The Kronecker product offers many properties some of these hsted below (for further 
details and proofs see Golub and Van Loan [9 , Van Loan [5S], A.N. Langville, W.J. 
Stewart 14J): 

1. Associativity. 

A(S>{B®C) = {A®B)®C) 

2. Distributivity. 

{A + B)®{C + D) ^ A®C + B®C + A®D + B®D 

3. Compatibility with ordinary matrix multiplication. 

AB®CD = {A® C){B ® D) 

4. Compatibility with ordinary matrix inversion. 

{A®B)-^ = A-'^ ®B-^ 

5. Compatibility with ordinary matrix transposition. 

{A (g) B)^ = A'^ (E)B'^ 



6. Trace factorization 



7. Norm factorization 



tr{A®B) = tr{A)tr{B) 



\A®B\\ = \\{A)\\\\{B)\\ 



8. Compatibility with Cholesky decomposition. 

Let A and B semi-definite positive matrices then: 

A®B^ [CaCD {CbCD = [Ca Cb){Ca Csf 

9. Special matrices. 

Let A and B be nonsingular, lower (upper) triangular, banded, symmetric, posi- 
tive definite, . . . , etc, then A® B preserves the property. 

10. Eigenvalue and Eigenvectors. 

Define two square matrices A and B, N x N and M x M, respectively. Suppose 
that Ai, . . . , Aat G ct{A), vi, . . . , vn and /^i, . . . , /j,m G ^(B), wi, . . . , wj^ are the 
eigenvalues and the correspondent eigenvectors of the two matrices respectively, 
where cr(.) denotes the spectre of the matrix. The Kronecker product. A® B, 
has eigenvectors Vi (g) wj and eigenvalues A^/ij. 

Summarizing, every eigenvalue of A g) _B arises as product of eigenvalues of A 
and i?, and every eigenvector as a Kronecker product between the corresponding 
eigenvectors. This last property still holds for singular value decomposition. 
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6 Generating Sample Path 



In discussing the simulation of a geometric brownian motion we should focus on the 
realization of a simple brownian motion at the sample time points of the grid. 

Because brownian motion has independent and normally distributed increments, 
simulating Wi{ti) is straightforward. 

Let ei, . . . , Eat be independent standard normal random variables and set Wi{to) — 0. 
Subsequent values can be generated as follow : 

W^{ti) = W^iU^i) + ^ti-ti^iei,l = l,...,N (33) 

For a brownian motion Xi{t) =BM(/ii, cji) given Xilto) set 

X,{ti) = X,{ti^i) + {tl - <(_i) + ^ti-ti^ia,€ul = 1, . . . , TV (34) 

For time-dependent parameters the recursion becomes (in the general situation the drift 
can be time-dependent too): 

X,{ti) = X,{ti^i) + ' (s) ds + ^JJ\f{^sei,l = 1,...,N (35) 

The methods (|33|) -(|35 |) are exact in the sense that the joint distribution of the random 
vector {Wi{ti), . . . , Wi{tN)) or {Xi{ti), . . . ,Xi{tN)) coincides with that of the original 
process at the times {^i, . . . , iAr}, but are subject to a discretization error. 

Nothing can be said about what happens between the time point of the grid. One 
might choose a linear interpolation to get intermediate values of the simulated process 
without obtaining a correct joint distribution. 

Applying the Euler scheme for the brownian motion with time-dependent drift and 
diffusion, 

X,{ti) = X,{ti^i) + ^i,{ti) {tl - ti^i) + ^ti-ti^ia,{ti)ei,l = 1,...,N (36) 

we introduce a dicretization error even at time points {^i, . . . ,ijv}, because the incre- 
ments will no longer have the correct mean and variance. 

The vector {Wi{ti), . . . ,Wi{tpf)) is a linear combination of the vector of the incre- 
ments {Wi{ti) — Wi{to), . . . , Wi{tN) — Wi{tN-i)) that is normally distributed. AU lin- 
ear combinations of normally distributed random vectors are still normally distributed. 

In general, let Y = CX be a A'^-dimensional random vector with multi-dimensional 
distribution N{^y, Sy) written as a A^x M linear transformation C of a M-dimensional 
random vector X with multi-dimensional distribution N{^x, Sx) then: 

Ey = CYsxC^. (37) 

This result provides an easy way to generate a vector of dependent normal random 
variables Y = CX ~ A^(/iy,Ey) from a set of independent ones X. Indeed, the 
dependence is completely taken into account by the covariance matrix: 

Sy = CC^ (38) 

The general problem consists of finding the linear transformation C, (for further details 
and proofs see Cufaro-Petroni 0). 

6.1 Cholesky Construction 

As far as the generation of a brownian motion is concerned, we note that method p3|) 
can be written as: 

/ W,{ti) \ / ^1 \ 

U ; , (39) 

V {tn) I ) 
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where Cr is the Cholesky matrix associated to the autocorrelation matrix of each 
brownian motion Wi (t) . 

Referring to the general problem the Cholesky decomposition simply faces the ques- 
tion of finding a matrix fulfilling equation (1381) among all lower triangular matrices. 

This is not a unique possibility, there are several other choices, but all of them must 
satisfy the general problem We will concentrate on two of them: the Principal 

Component Analysis (PCA) proposed by Acworth, Broadie, and Glasscrman (1998) 
[1] and a Kronecker Product Approximation that we introduce as a different and new 
approach in Section 5.4. 

We apply the Cholesky decomposition method in order to draw the random vector 
e with distribution iV(0, Emat)- 

In case of constant volatilities we showed that T,mn — R<E)'S,. We can exploit the 
Kronecker product compatibility with Cholesky decomposition to get: 

C^MN = Cr (g> Cs (40) 



where Cr is given by equation ([29|) . By means of the Kronecker product we can reduce 
the computational effort by splitting the analysis of an AfA^ x MN matrix into the 
analysis of two smaller M x M and N x N matrices . 

When time-dependent volatilities are considered we cannot exploit the properties of 
the Kronecker product. Sa/jv can be partitioned into block matrices S(ii), . . . , S(i7v) 
that are not constant anymore and depend on the point of the time grid. 

Provided this time-dependent feature, all the information carried out by Smat hinges 
in N smaller M x M matrices. These latter matrices depend on the particular time- 
dependent functions that determine the evolution of the volatilities and on the constant 
correlation among the assets returns (the analysis can be applied to time-dependent 
instantaneous correlations). 

In the following we present a faster than the standard Cholesky decomposition 
algorithm that focuses on particular form of the covariance matrix Smat- 

In the time-dependent volatility case the global covariance matrix 'Smn satisfies the 
boomerang shape property as R as well as their Cholesky matrices. We consider this 
feature with respect to the partitioned matrix notation. 

It is possible to develop all the calculations storing A^ block matrices, (E(ti), . . . , E(tjv)), 
in a tri-linear tensor (Stot)ifci. For any fixed / the block (Stot)^^.; coincides with S(i;-). 
Consequently we perform the ad hoc Cholesky decomposition suited for partitioned 
boomerang shaped matrices. 

Using the partitioned matrix notation, the Cholesky algorithm develops according 
the following steps: 









^BR 



Ctl 





Cbl 


Cbr 



Ctl 


^BL 





Cbr , 



The block matrices with index TL (Top-Left) are M x A/, those ones with BL (Bottom- 
Left) are (A^- l)Af x Af, those ones with 5ir:(Bottom- Right) are {N -\)M x {N~\)M. 

1. Decompose the Top- left block. 

Sr^ /^T r~i /^T 

TL — '-^TL'-^TL — 

2. Decompose the Bottom- left block. 

S_BL — CblCtl 

In particular exploiting the boomerang shape property we should have: 

/ \ / CtlCT 



^BL 



TL 



\ ^TL / \ CtlC^j^ 

Due to the boomerang shape structure of the global covariance matrix, this second 
step can be avoided, because it consists of repeating the first step. 
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3. The Cholesky decomposition is iterated to Bottom-Right block. 



S_Bi?, — CbrCqr + CblCq]^ 

The last term on the right hand side of the previous equation is known, because 
it has been calculated in step 1. 

We let ^Update define a {N— 1)M x (N — 1)M matrix by the following expression: 



^Update — ^BR — CblC^i^ — LjBR'^BR 

we can conclude that after decomposing Tjjjpdate and getting we have the 
complete picture of the global Cholesky matrix. 

This last step can be specified in greater detail referring to the boomerang shape 
feature of T^update- 



CbrCT 



^Update 



nt2) 


^TR. 




^BR. 



I CT 



- G 



TL 



.C 



TL 



TL 



where Ebl and Es/j, are {N - 1)M x M and {N - 1)M x (TV - 1)M matrices. 
After all the calculation we obtain: 



-'Update 



2(^2) — CtlC^l 


TR 


BL 


BR 



CbrCbr 





TR 


BL 


BR 



where TR, BL and BR are partitioned boomerang shaped matrices. C2 represents 
the M X M Top-Left block of Cbr, while E(ii) = CtlC^l = CiCf 

The algorithm can be implemented running a loop of N iterations. 

The first iteration consists of realizing the Cholesky decomposition of step 1 de- 
scribed above. 

The generic iteration i consists in subtracting the Top-left Block of the i — l updated 
matrix to all the remaining N — i blocks (their dimension is M x M) of the tri-linear 
tensor (Etot)ifci and that calculate the calculate the Cholesky decomposition. 

This algorithm returns N block matrices, whose dimension is M x M, that are 
stored in tri-linear tensor, {Ctot)ikj that represents the global Cholesky matrix. 



6.2 Principal Component Analysis 

A more efficient approach for the path generation is based on the Principal Component 
Analysis (PC A). 

is a symmetric matrix and can be diagonalized as 

Sy = EAE^ = {EA^/^){EA^/^f . (41) 

For this method, the linear transformation C solving equation ([55]) is defined as Eh^^"^ . 
A is the diagonal matrix of all the positive eigenvalues of sorted in decreasing order 
and E is the orthogonal matrix {EE^ ~ 1) of all the correspondent eigenvectors. 

The matrix Eh^l"^ has no particular structure and generally does not provide com- 
putational advantage with respect to the Cholesky decomposition. 

This transformation can be interpreted as a sort of rotation of the random vec- 
tor whose covariance matrix is Sy; in the new frame of reference it has independent 
components whose variances are the elements on the diagonal of A. 

The higher efficiency of this method is due to the statistical interpretation of the 
eigenvalues and eigenvectors (see Glasserman jl]). 

Suppose we want to generate Y ^ A^(0, Ey) from a vector e ~ N{0,I), we know 
that the random vector can be set as: 

d 

Y = ^ CkCfc 

fe=i 
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where Ck is the k-th column of C. 

Assume has full rank d, then it is non singular and invertible and the factors 
efe are themselves linear combination of Yfe. In the special case C = EA^^^, ek is 
proportional to ek • Y. 

The factors ek constructed in the previous way are optimal in a precise statistical 
sense. 

Suppose we want to find the best singled- factor approximation of Y, that is to find 
the best linear approximation that best captures the variability of the components of Y. 
The optimization problem consists in maximizing the variance of w • Y with constraint 
of the form w • w = 1 : 

max w • Eyw (42) 

w-w=l 

If we sort the eigenvalues of Sy in decreasing order then the optimization problem is 
solved by ei. More generally the best fc-factors approximation of Y leads to factors 
proportional to ei • Y, . . . , ek ■ Y with ei • e^ = Sim, with: 

efc = -^ek • Y. (43) 
V Afe 

This representation can be recasted as the minimization of the mean squared error: 



MSE = E 



1 



(44) 



where we are looking for the best fc-factors mean square approximation of X. This 
formulation gives the same results. 

In the statistic literature the linear combination ek • Y is called principal component 
of Y. The amount of variance explained by the first k principal components is the ratio: 

where d is the rank of Ey. 

We can apply PCA to generate a onc-dimcnsional brownian motion BM(0, R) cal- 
culating the eigenvectors and eigenvalues of the sampled auto-covariance matrix R and 
then rearranging them in decreasing order. The magnitude of the eigenvalues of this 
matrix drops off rapidly. For instance it is possible to verify that in the case of a brown- 
ian motion with 32 time steps the amount of variance explained by the first five factors 
is 81% while it exceeds 99% at fc = 16. 

This result is fundamental in identifying the effective dimension of the integration 
problem. PCA helps Monte Carlo estimation procedures based on the generation of 
brownian motion where we should identify the effective dimension of the problem. With 
this choice we can identify the most important factors in a precise statistical framework 
by fixing a value p in the determining the effective dimension, (for instance p = 99%). 

This statistical ranking of the normal factors cannot be implemented by Cholesky 
decomposition that we expect will return unbiased Monte Carlo estimations but higher 
RMSEs. 

As far as the multi-dimensional brownian motion is concerned, we start with the 
constant volatility case. We have already shown in section 4 that the covariance matrix 
J^MAT of the multi-dimensional brownian motion BM(0, Yimn) can be written as i?® S. 

Property 10 of the Kronecker product permits to improve the speed of the compu- 
tation of the eigenvalues and eigenvectors of Sa/at. It reduces this calculation into the 
computation of the eigenvalues and vectors of the two smaller matrices R and S. 

Coupling the use of the Kronecker product analysis with the ANOVA definition of 
effective dimension we can implement a fast and efficient Monte Carlo estimation in 
order to price exotic multi-dimensional path-dependent options. 

Empirical evidence in finance shows that effective dimension is often lower than the 
problem dimension d, (see Caflisch, Morokoff, and Owen 4\ for a general discussion). 
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We focus our analysis on Asian options pricing after formulating the pricing problem 
an integral. As we presented in section 3 ANOVA is used as to provide a representation 
of the integrand as a sum of orthogonal functions. If each of these orthogonal functions 
depends only on a distinct subset of the coordinates, the integrand can be written as 
a sum of integrals of functions of lower dimension. The complexity of the computation 
of the integral has been reduced with respect to the integral dimension. In pricing 
Asian options we are not able to reduce the dimension of the original integrand by this 
approach, because we cannot exactly find a set of orthogonal functions. What we can 
propose is an approximation based on the PCA construction. In our finance problems 
we achieved a representation involving matrices, describing the dependence between the 
different variables, as arguments of the exponential function <?(.). Our approximation 
consists in a direct application of ANOVA and effective dimension calculation to the 
random vector Z. This is equivalent to the Taylor expansion up to the first order of the 
exponential function g(.) that leads to the following definition of effective dimension, 
dx, of the problem (in truncation sense): 

dT 

^ Arf < tr{A)p (46) 

d=l 

where Xd € ct(Sa./jv). The level p is arbitrary; we chose p — 99%. 
6.3 The Kronecker Product Approximation 

The time-dependent volatilities market has a covariance matrix Smtv with time-dependent 
blocks. Generally, it has not a particular expression because it depends on the volatility 
functions and the instantaneous correlation. The covariance matrix of the asset returns 
is not anymore constant so that cannot be written as a Kronecker product. 

We have shown that a fast Cholesky decomposition algorithm can be ran but it does 
not take any ANOVA and effective dimension consideration, while the PCA approach 
is still applicable but we cannot reduce the computational burden using the properties 
offered by the Kronecker product. 

In the constant volatility case the special structure of T,mn makes possible to com- 
pute all the eigenvalues and eigenvectors with M^+N^ operations, written 0{M^ +N^), 
instead of O ((MiV)'^) for a general MN x MN square matrix. 

The market under consideration has the multi-dimensional brownian motion as 
unique source of risk. Its generation procedure is independent of the constant or time- 
dependent volatilities because its autocovariance matrix R is not influenced by these 
market features. 

Based on these considerations our proposition is to find a constant covariance ma- 
trix among the assets, K, in order to approximate, in an appropriate sense, the global 
covariance matrix as a Kronecker product of R and K. Our hypothesis is that the 

effective dimension of the problem should not dramatically change after this transfor- 
mation with an advantage from the computational point of view. We develop the PCA 
decomposition of the approximating matrix assuming that the principal components are 
not so different from those of the original random vector. This approximation would 
lead to a different multi-dimensional path because Ri^ K is not the covariance matrix 
of the original process. The global and true path is reobtained using the Cholesky 
factorization. 

In the following we illustrate the proposed procedure that we label KPA. 

The general problem consists of finding two matrices B S K'"!^*"! and C G K™2><"2 
that minimize the Frobenius norm. All calculations and proofs can be found in Pitsianis, 
Van Loan f22^ and Van Loan [25] : 

^AiB,C) ^\\ A- B^C \\^ (47) 

where A G K™^" is an assigned matrix with m = niim2 and n = nin2- 

The main idea is to look for a rearrange matrix 'R-{A) such that equation (|T7)) can be 

rewritten as $^(5, C) =\\ TKyA) - vec{B) (g) vec{CY f • 
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Definition 5 (The vec operation). The vec operation transforms a matrix X e R*^'^ 
into a column vector vec{X) G M''^^"'^ by 'stacking' the columns: 



' 021 022 ' 



021 

ai2 

V 022 J 

As far as our approximation is concerned the general problem is simplified. Indeed, 
the new problem consists of finding only one matrix K minimizing the Frobenius norm: 

<i>{K)^\\j:MN~R<E)Kf (48) 

The approach is equivalent to a Least Square problem in the Kik. 

The elements Kik are given by the formula below (for a complete proof see Pitsianis, 
Van Loan [12] p. 8): 

tr(TZ{'SMN)ikR) 

K^k = ^ (49) 

tr RRT^^ 



where TZ{'SMN))ik is a N x N matrix. For any i and k ranging from 1 to M, TZ{YiMN))ik 
is obtained by sampling Y.mn with M as sampling step. 

By its definition, it can be noticed that for any i and k TZ{T,MN))ik is a boomerang 
shaped block matrix. 

By direct computations and relying on the particular form of i?, the denominator 
of the equation ([35]) is: 

N 

tr[RR^^ ^tr[R^^ = ^ {2{N - j) + \^t] (50) 

Moreover, given two general N x N boomerang shaped matrices A and B the trace 
of their product is: 

N 

tr{A^B) = tr{AB) = ^ (2(iV - j) + l)a,,5,, (51) 

Ojj and bjj are the only significant value to store. 

The considerations above permit to evaluate X in a fast and efficient way without 
high computational efforts. 

As already mentioned, if we would use the ANOVA-PCA procedure to i? (g) X we 
would not get the required path. Let E and A be the eigenvectors and eigenvalue 
matrices associated to R® K^ if we would consider EA}/'^ as a generating matrix we 
would generate a path whose global covariance matrix 'is R® K and not S a/at . 

In order to tackle to the original problem the Cholesky decomposition is used. In 
fact given two N dimensional random vectors Zi and Z2 with covariance matrices Si 
and S2 respectively, we can always write: 

HZ <^^' 

where Ci and C2 are the Cholesky matrices of Si and E2, respectively and e is a vector 
of independent random variables. At the same time we can generate Z2 by PCA: 

Z2 = E2Ky\ (53) 

where E2 and A2 comes from the complete PCA of S2. 
Combining the above equalities we have: 

Zi = CiC^^E2Ky^e (54) 
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It is possible to generate a random path Zi applying the PC A to Z2 and than turn- 
ing back to the original problem. Our fundamental assumption is that the effective 
dimension of our problem remains almost unchanged and, in the estimation procedure, 
we apply almost the same statistical importance to the original principal components 
giving an advantage from the computational point of view. 

Focusing this result to the problem under study, we let Si — Yimn and E2 — R®K 
so that equation (f54|) becomes: 



-^2^i- ' e (55) 

Cemn) ^r^ and are the Cholesky matrices of Smat, R and K, respectively. In 
the derivation of the previous equation we exploit several properties of the Kronecker 
product. 

We again stress the fact that in the case of time-dependent volatilities we analyze the 
effective dimensions of the integral problem after the Kronecker product approximation. 
Generally this second approximation would return a higher effective dimension with 
respect to the normal case where only a linear approximation is considered. Furthermore 
our method generates the correct required path as proved. 

In order to obtain a fast and efficient algorithm for the path generation we develop 
all the calculations: 



1. Cn^ ® Cj^^ is a sparse bi-diagonal partitioned matrix: 



CI 



K 



%/tr 








\ 



2. C'sJ^J„Cp^ (E) Cj^^ is lower triangular partitioned matrix. 



c 



K 



C^3 



CI 



C^i for i = I, . . . , N indicates the i-th block matrix of the tri-linear tensor {Ctot)ikj- 
Ai = ti — where = is understood. 

Only {Ctotjikj and the sequence { Ai}i=i^...^7v need to store all the information em- 
bedded in Csm„C^^ ® C^. 

The total generating matrix C-£,^ff^C^ ® C^^£'2A^/^ can be computed quickly by ma- 
trix product with partitioned matrices. 



7 Solution Methodology 

We aim to provide an efhcient technique that improves the precision of the general 
Monte Carlo method to exotic derivative contracts and in particular Asian options. 
According to equation (|23p the actual problem consists of generating a sample of uni- 
form random draws to uniformly cover the whole hypercube [0, 1] . In the following 
subsections we introduce different ways to generate random numbers that uniformly 
cover the hypercube [0, l]''. 
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7.1 Stratification and Latin Hypercube Sampling 

Stratified sampling is a variance reduction method for Monte Carlo estimates. It 
amounts to partitioning the hypercube T) = [0, 1)'' into H disjoint strata {h — 
1, . . . ,H), i.e., V = UiLi T^h where = for all j ^ k, then estimating the 

integral over each set, and finally summing up these numbers (see Boyle, Broadie and 
Glasserman [3] for more on this issue). Specifically, mutually independent uniform 
samples Xi, . . . , are simulated within a stratum T^h , and the resulting integrals are 
combined. The resulting stratified sampling estimator is unbiased. Indeed: 



71=1 i=l 
H 

h=l 
H 



J2 f n^)dx^ I. 
h=i 



where \'Dh\ denotes the volume of stratum Dh. Moreover, this estimator displays a 
lower variance compared to a crude Monte Carlo estimation, i.e., 



Var 



I strat 



< 



Stratified sampling transforms each uniformly distributed sequence Uj = {Uij, . . . , Udj) 
in V into a new sequence Vj = (Vij, . . . , Vdj) according to the rule 



\Jj + (ii, ...,id) 



= 



, n,ifc = 0, . 



l,k^ 1,. 



where (ii, . . . ,id) is a deterministic permutation of the integers 1 through d. This 
procedure ensures that one Vj lies in each of the n'^ hypercubes defined by the strati- 
fication. Latin Hypercube Sampling (LHS) can be seen as a way of randomly sampling 
n points of a stratified sampling while preserving the regularity from stratification (see, 
for instance, Glasserman [5]). Let tti, . . . jiTd be independent random permutations of 
the first n positive integers, each of them uniformly distributed over the n\ possible 
permutations. Set 



Ujk + TTfc (j) - 1 



j = 1, . . . ,n,fc = 1, 



(56) 



where tt^ (j) represents the j-th component of the permutation for the fc-th coordi- 
nate. Randomization ensures that each vector Tj is uniformly distributed over the d 
dimensional hypercube. Moreover, all coordinates are perfectly stratified since there is 
exactly one sample point in each hypercube of volume 1/n. For d = 2, there is only 
one point in the horizontal or vertical stripes of surface 1/n (see Figure 2). The base 
and the height are 1/n and 1, respectively. For d > 2 it works in the same way. It can 
be proven that for all n > 2, d > 1 and squared integrable functions /, the error for the 
estimation with the Latin Hypercube Sampling is smaller or equal to the error for the 
crude Monte Carlo (see Koehler and Owen p^): 



Var 



Ilhs 



< 



1 



(57) 



Figure 2 shows the distribution of 32 points generated with the LHS method. For the 
LHS method we notice that there is only 1 point (dotted points in Figure 2) in each 
vertical or horizontal stripe whose base is 1 and height is 1/32: it means that there is 
only a vertical and horizontal stratification. 
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CRUDE MC 



« 



Figure 1: The panel shows 32 points drawn with standard pseudorandom generators 




Figure 2: The panel shows 32 points generated with LHS 



19 



7.2 Low-Discrepancy Sequences 

As previously mentioned, the standard MC method is based on a completely random 
sampling of the hypercube [0, 1)'' and its precision can be improved using stratification 
or Latin Hypercube sampling. These two methods ensure that there is only one point 
in each smaller hypercube fixed by the stratification as illustrated in Figure 2. At 
the same time, these techniques provide nothing more than the generation of uniform 
random variables in smaller sets. 

A completely different way to approach the sampling problem is to build-up a de- 
terministic sequence of points that uniformly covers the hypercube [0, 1)'^ and to run 
the estimation using this sequence. Obviously, there is no statistical quantity that 
may represent the uncertainty since the estimation always gives the same results. The 
Monte Carlo method implemented with the use of low-discrepancy sequences is called 
Quasi-Monte Carlo (QMC). 

The mathematics involved in generating a low-discrepancy sequence is complex and 
requires the knowledge of the number theory. In the following, only an overview of the 
fundamental results and properties is presented (see Niederreiter [17j for more on this 
issue) . 

We define the quantity D* = D* (Pi, . . . , Pn) as the star discrepancy. It is a measure 
of the uniformity of the sequence {-Pn}„gN* G [0, 1)'' and it must be stressed that it 
is an analytical quantity and not a statistical one. For example, if we consider the 
uniform distribution in the hypercube [0, 1)'^ , the probability of being in a subset of 
the hypercube is given by the volume of the subset. The discrepancy measures how 
the pseudo-random sequence is far from the idealized uniform case, i.e. it is a measure, 
with respect to the L2 norm for instance, of the inhomogeneity of the pseudo-random 
sequence. 

Definition 6 (Low-Discrepancy Sequencies). A sequence {Pn}„gN* is called low- 
discrepancy sequence if: 

D:(Pi,...,P„) = o(^^i^^. (58) 

i.e. if its star discrepancy decreases as {\nn)'^ /n. 

The following inequality, attributed to Koksma and Hlawka, provides an upper 
bound to the estimation error of the unknown integral with the QMC method in terms 
of the star discrepancy: 

\I-i\<D:VHK{.f). (59) 

Vhk if) is the variation in the sense of Hardy and Krause. Consequently, if / has a 
finite variation and n is large enough, the QMC approach gives an error smaller than the 
error obtained by the crude MC method for low dimensions d. However, the problem 
is difficult owing to the complexity if estimating the Hardy-Krause variation, which 
depends on the particular integrand function. 

In the following sections we briefly present digital nets and the well-known Sobol' 
sequence that is the most frequently used low discrepancy sequence to run Quasi-Monte 
Carlo simulations in finance. 

7.3 Digital Nets 

Digital nets or sequences are obtained by the number theory and owe their name to the 
fact that their properties can be recognized by their digital 6-ary expansion in base b. 
Many digital nets exist; the ones most often used and considered most efficient are the 
Sobol' and the Niederreiter-Xing sequences. 

The first and simplest digital sequence with d = 1 is due to Van der Corput and is 
called the radical inverse sequence. Given an integer 6 > 2, any non-negative number 
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N 


n base 2 


(1)2 {n) base 2 


4>2 (n) 





000. 


0.000 


0.000 


1 


001. 


0.100 


0.500 


2 


010. 


0.010 


0.250 


3 


oil. 


0.110 


0.750 


4 


100. 


0.001 


0.125 


5 


101. 


0.101 


0.625 


6 


110. 


0.011 


0.375 


7 


111. 


0.111 


0.875 



Table 1: Van der Corput sequence. 



n can be written in base b as: 

oo 

n = ^nfe6'=-i. (60) 

fc=i 

The base b radical inverse function (n) is defined as: 

oo 

cj)t{n) = J2^kb-'' e[0,l), (61) 
fe=i 

where nfc G {0, 1, . . . , 6 — 1} (Galois set). 

By varying n the Van der Corput sequence is constructed. Table [T] illustrates the 
first seven Van der Corput points for 6 = 2. Consecutive integers alternate between odd 
and even; these points alternate between values in [0, 1/2) and [1/2, 1). The peculiarity 
of this net is that any consecutive 6™ points from the radical inverse sequence in base 
b are stratified with respect to b™ congruent intervals of length 6^™. This means that 
in each interval of length b~™ there is only one point. 

Table [T] shows an important property that is exploited in order to generate digital 
nets, because a computing machine can represent each number with a given preci- 
sion, referred to as "machine epsilon". Let z = O.Z1Z2 ■ ■ ■ (base b) G [0,1) , define 
'4'(z) = (zi,Z2, ■ ■ • ) the vector of the its digits, and truncate its digital expansion at 
the maximum allowed digit w: z = J2k=i ^kb^*' ■ Let n = \b'^z\ = X)r=i ^hb'^^^ £ N* , 
where [x] denotes the greatest integer less than or equal to x. It can be easily proven 
that: 



nh 



<^w-h+l 



(z) V/l = 1, 



This means that the finite sequences {nh}f^^^i and {zk}i^^^i have the same 
elements in opposite order. For example, in the table [T] we allow only 3 digits; in order 
to find the digits of 02 (1) = 0, 5 we consider (/)2 (1) 2^ = 4 = 0?ii + ?i20 + n^l. The 
digits of 4>2 (1) are then (1, 0, 0) as shown in the table [TJ 

The peculiarity of the Van der Corput sequence is largely required in high dimen- 
sions, where the contiguous intervals are replaced by multi-dimensional sets called b-adic 
boxes. 

Definition 7 (b-iadic Box). Let b > 2, kj, Ij with < Ij < b'^^ be all integer numbers. 
The following set is called b-iadic box: 



d r 



n 



L L + 1 



(62) 



where the product represents the Cartesian product. 

Definition 8 ((t,m,d) Nets). Let t < m be a non-negative integer. A finite set of 
points from [0, 1)'' is called {t,m, d)-net if every b-adic box of volume 5"™+* (bigger 
than b^"^ ) contains exactly 6* points. 
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d 


P 


M 


Principal polynomial 


q 


1 


[1] 


[1] 


1 





2 


[1 1] 


[1] 


x + 1 


1 


3 


[1 1 1] 


[1 1] 


+ X + 1 


2 


4 


[10 11] 


[1 3 7] 


+ X + 1 


3 


5 


[110 1] 


[1 1 5] 


x^ + x^ + 1 


3 


6 


[10 11] 


[13 11] 


X^ + X + 1 


4 


7 


[110 1] 


[113 7] 


x"^ + X'^ + 1 


4 


8 


[10 10 1] 


[1 3 3 9 9] 


x^ + x^ + 1 


5 


9 


[1110 11] 


[1 3 7 13 3] 


X^ + X^ + X^ + X + 1 


5 


10 


[10 1111] 


[1 1 5 11 27] 


X^ + X^ + X^ + X + 1 


5 



Table 2: Initial values satisfying Sobol property A up to dimension 10. By convention, the 
recurrence relation for the 0-degree polynomial is M^, = 1 

This means that cells that "should have" 6* points do have 6* points. However, 
considering the smaller portion of volume b~"^, it is not guaranteed that there is just 
one point. 

A famous result of the theory of digital nets is that the integration over a [t, m, d) 
net can attain an accuracy of the order of O ^In''"^ (n) /nj while, restricting to {t,d) 

sequences, it raises slightly to O (lii^ (n) /nj (see Niederreiter The above results 

are true only for functions with bounded variation in the sense of Hardy-Krause. 

7.4 The Sobol' Sequence 

The Sobol' sequence is the first d dimensional digital sequence, (6 = 2), ever realized. 
Its definition is complex and is covered only briefly in the following. 

Definition 9 (The Sobol' Sequence). Let {"-fcjfcgpj. he the set of the h-ary expansion 
in base b ~ 2 of any integer n; the n-th element Sn of the SohoV sequence is defined as: 

+ C30 

5„ = ^(nfcFfcmod2) 2-^ (63) 
fc=i 

where 14 G [0, l)'^ are called direction numbers. In practice, the maximum number of 
digits, w, must be given. In Sobol's original method the «-th number of the sequence 
Sij, i e N, j e {1, . . . , d}, is generated by XORing (bitwise exclusive OR) together the 
set of Vkj satisfying the criterion on k : the fc-th bit of i is nonzero. Antonov and Saleev 
derived a faster algorithm by using the Grey code. Dropping the index j for simplicity, 
the new method allows us to compute the {i + l)-th Sobol' number from the i-th by 
XORing it with a single Vk, namely with fc, the position of the rightmost zero bit in i 
(see, for instance. Press, Teukolsky, Vetterling and Flanncry 23] ). Each different Sobol' 
sequence is based on a different primitive polynomial over the integers modulo 2, or 
in other words, a polynomial whose coefficients are either or 1. Suppose P is such a 
polynomial of degree q: 

P = + aix'"^ +a22:«"2 H ha^-ix + l. (64) 

Define a sequence of integers Mfe, by the qth term recurrence relation: 

Mk = 2aiMfc_i e 2^a2Mk-2 ® • • • ® ® (2'?Mfc_^ Mk-q) ■ (65) 

Here © denotes the XOR operation. The starting values for the recurrence are Mi, . . . , Mg 
that are odd integers chosen arbitrarily and less than 2, . . . , 2^, respectively. The direc- 
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Sobol 




Figure 3: The panel shows the first 32 points of the 2-dimensional Sobol' sequence 



tional numbers Vk are given by: 

Vk^^ k^l,...,w. (66) 

Table [5] shows the first ten primitive polynomials and the starting values used to gen- 
erate the direction numbers for the 10 dimensional Sobol' sequence. 

7.5 Scrambling Techniques 

Digital nets are deterministic sequences. Their properties ensure good distribution in 
the hypercube [0, 1)'^, enabling precise sampling of all random variables, even if they are 
very skewed. The main problem is the computation of the error in the estimation, since 
it is difficult to compute and depends on the chosen integrand function. To review, 
the crude MC provides an estimation with low convergence independent of d and the 
possibility to statistically evaluate the RMSE. On the other hand, the QMC method 
gives a higher convergence, but there is no way to statistically calculate the error. 

In order to estimate a statical measure of the error of the Quasi-Monte Carlo method 
we need to randomize a (t, m, c?)-net and try to obtain a new version of points such that 
it still is a {t,m,d)-net and has uniform distribution in [0, 1)''. 

This randomizing procedure is called scrambling. The scrambling technique per- 
mutes the digits of the digital sequence and returns a new sequence that has both the 
properties described above. 

The scrambling technique we use is called Faure-Tezuka Scrambling (for a precise 
description see Owen [21], Hong and Hickernell flOj). 

For any z E [0, 1) we define ^{z) as the oo x 1 vector of the digits of z. 

Now let Li, . . . be nonsingular lower triangular cx) x oo matrices and let ei, . . . , 
be oo X 1 vectors. Only the diagonal elements of Li, . . . are chosen randomly and 
uniformly in = {1, . . . , 6} , while the other elements are chosen in Zt, = {0, 1, . . . , 6}. 
Y, the Faure-Tezuka scrambling version of X, is defined as: 

* (y^j ) = (ij * (xy ) + ) modb (67) 
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Figure 4: The panel shows the first 32 point of the Sobol sequence compared to their 
Faure-Tezuka scrambled version 



All operations take place in the finite field Zi,. Owen proved that, with his scrambling, 
it is possible to obtain (see Owen [18j): 





6* 






< — 






n 


[h-lj 



(68) 



for any twice integrable function in [0, 1) . These results state that for low dimension 
d , the randomized QMC (RQMC) provides a better estimation with respect to Monte 
Carlo, at least for large n. 



8 Implementation and Algorithm 

We illustrate the simulation procedure to compute the arithmetic Asian option price. 
The purpose of our analysis is to characterize the efficiency of Monte Carlo methods 
based on the path generation techniques and the uniform points used for the evaluation 
of the integral We consider separately the constant volatility and time-dependent 
volatility markets. 

It must be stressed that Quasi-Monte Carlo estimations are dramatically influenced 
by the problem dimension, because the rate of convergence depends on the problem 
dimension d, as it can be seen in equations (|58p and (|68p . Many studies and experiments 
suggest that Quasi-Monte Carlo methods can only be used for problem dimensions up 
to 20 (see Boyle, Broadie and Glasserman [2] for more on this issue). This condition 
translates into a relationship between the number M of underlying assets and the 
number N of monitoring times: M x A'^ < 20. When this condition is not satisfied 
anymore we use the Latin Supercube method that we describe hereafter. 

8.1 Latin Supercube Sampling 

The scrambling procedure allows the statistical estimation of the RMSE as the crude 
MC does with the order of convergence that depends on the the dimension d. For high 
d the fast convergence of the RQMC is lost, there is no benefit to use it compared to the 
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simple MC. Generally in finance the dimension is high even using dimension reduction 
techniques like ANOVA-PCA decomposition. 

Owen [5D] has proposed a method to extend the convenience of applicability of RQMC 
for high dimensions. This method is called Latin Supercube Sampling, (LSS), owing to 
its similarity to the LHS. The random permutation is now applied to a set of subse- 
quence of the original one with some statistical sense. 

Let Y = {yi, . . . , yb""} be the digital sequence of the simulation variables, and 6™ = N. 
Dividing it into k nonempty and disjoint subsets Y = IJ^^j^ Yr and letting Sr = dirriYj- 
we have X]r=i — ^- practice, each point of the sequence can be represented as 
yi ~ (Xi 1 ■ • • iXi)^ where Xi ^ [0, these points Xi £^re ordinarily points of an Sr- 
dimensional RQMC method. 

For r — l,...,fc let 7rr(i) be an independent uniform and random permutation of 
{!,..., N} than a Latin Supercube sample is obtained by taking: 

yi = (xii«,---,x',w) (69) 

It means that the first si colums in the LSS are obtained by randomly permuting the run 
order of the RQMC points xh ■ ■ ■ iXh the next S2 columns come from an independent 
permutation of the run order of xf so on. 

The convenient way to divide the original set might be arranging them in statistically 
orthogonal sets using the ANOVA-PCA decomposition. 

In practice in financial simulation with d Brownian motions, it may make sense to 
select 5 principal components of each path, to apply an RQMC method to each of them 
with LSS and then pad them out the other variables with LHS. In fact, about 95% of the 
total variance of the Brownian motion is explained by these components. Alternatively, 
it may be better to group the first k principal component, then the second, and so on. 

However, all these results are weak and only the practical test can give an answer 
to which sequence and scrambling should be used. 

8.2 Key Steps of the Simulation Procedure 

As a first scenario we run simulations using the Cholcsky and the PCA decomposition 
procedures for the constant volatility case. As a second scenario we test the efficiency 
of the proposed Kronecker product approximation by comparing the its results with 
those obtained with the PCA decomposition. 

As a random number generator we use three configurations: standard, LHS and 
Faure-Tezuka scrambled version of the Sobol' sequence. 

The test for constant volatility consists of three main steps: 

1. Random number generation by standard MC, LHS or RQMC. 

2. Path generation with Cholesky and PCA decompositions. 

3. Monte Carlo estimation. 

For the time-dependent volatility case the three steps are: 

1. Random number generation by RQMC. 

2. Path generation with PCA decomposition and KPA. 

3. Monte Carlo estimation. 

The first step of both cases is realized by using the correspondent random generator 
of uniform random variables. In order to extract normal random variables we rely on 
the inverse transform method that require the numerical inversion of the cumulative 
function of the standard normal. This numerical procedure may destroy the better 
stratification and the uniformity introduced by LHS and especially by low-discrepancy 
sequences. We use the Moro's algorithm that is more precise than the standard one 
due to Beasley and Springer. It provides a better accuracy on the tails of the inverse 
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Si(0) 


= 100 


K 


= 100 


r 


= 2% 


T 


= 1 


0"! 


= 30% 




= 40% 


Pij 


= and 40% fori,j = l,2. 



Table 3: Input Parameters Used in the First Simulation 

normal where we require that the LHS and Sobol sequences must reveal their higher 
precision, see Moro [16j and Glasserman ^ for more on the topic. 

For constant volatilities the second step can be implemented by the following algo- 
rithm: 

1. Define the parameters of the simulation. 

2. Define the drift as in equation ([201) . 

3. Create the N x N correlation matrix {R)^ — {ti /\ tk) ;l, k — 1, . . . ,N. 

4. Define the correlation matrix E among the M asset returns. 

5. Perform either a PCA or the Cholesky decomposition on the global correlation 
matrix T,mn- This matrix is built up by repeating the constant block of correla- 
tion S at all the times of observation. 

For time-depending volatilities we define the drift as equation ([21]) . while the last op- 
eration consists of performing the PCA decomposition and the KPA. 

Stratification introduces a correlation among random drawings so that the hypoth- 
esis of the Central Limit theorem is not satisfied and we cannot compute the RMSE 
straightforward. We rely on the batch method that consists of repeating Nb simulations 
for B times (batches). We assume that each of the B batches eliminates the correlation 
and the results form a sequence of B independent random variables. We compute the 
average Asian price for each batch; the RMSE becomes: 

where {a (0)^ , . . . ,a (0)^) is a sample of the average present values of the Asian option 
generated in each batch. 

9 Numerical Experiments 

We perform a test of all the valuation procedures described in the previous section. 
We specify our investigation into constant and time-dependent volatilities cases while 
our experiments involve standard Monte Carlo, the Latin Hypercube Sampling and 
Randomized Quasi Monte Carlo with the Faure-Tezuka scrambled version of the Sobol' 
sequences. 

9.1 Constant Volatility: Results 

As for a first pricing experiment we consider an at-the-money arithmetic Asian option 
with strike price K = 100, written on a basket of M = 2 underlying assets, expiring at 
T = 1 year and sampled N = 5 times during its lifetime. 

All results are obtained by using S — 8192 drawings and 10 replications. Table [3] 
reports the input parameters for our test. The nominal dimension of the problem is 
M X N ^ 10 that is equal to the number of rows and columns of the global correlation 



26 





Standard MC 


LHS 


RQMC 


PCA 
Cholesky 


7.195 (0.016) 
7.242 (0.047) 


7.157 (0.013) 
7.179 (0.022) 


7.1696 (0.0017) 
7.1689 (0.0071) 



Table 4: Uncorrelation Case. Estimated Prices and Standard Errors. 





Standard MC 


LHS 


RQMC 


PCA 
Cholesky 


8.291 (0.053) 
8.374 (0.055) 


8.2868 (0.0073) 
8.293 (0.026) 


8.2831 (0.0016) 
8.2807 (0.0064) 



Table 5: Correlation Case. Estimated Prices and Standard Errors. 



matrix T,mn- Paths are simulated by using both PCA and the Cholesky decomposition 
as in Dahl and Benth [S] and [7]. 

Table |4] and Table [5] show the results for the positive correlation and uncorrelated 
cases, respectively. Simulated prices of the Asian basket options are in statistical ac- 
cordance, while the estimated RMSEs depend on the sampling strategy adopted. The 
rate of convergence of the RQMC estimation is higher than the other two methods. 
In particular it is ten times higher than the standard Monte Carlo method that would 
return the same accuracy with 100 x S drawings. 

We observe that the PCA generation provides a better estimation both for LHS 
and RQMC, because these ones are more sensitive to the effective dimension, while 
PCA causes no distinction for the standard MC. The effect is more pronounced for 
the correlation case where the more complex structure of the global correlation matrix 
TiMN influences the estimation procedure. 

As from a financial perspective, it is normal to find a higher price in the positive 
correlation case than in the uncorrelated one. 

Moreover, we develop our analysis by investigating a very high-dimensional pricing 
problem. A basket of M = 10 underlying assets is considered with N — 250 sampling 
time points, the nominal dimension is d — 2500. 

We run our simulation with the same parameters used by Imai and Tan [llj and 
use the LSS for the high-dimensional QMC estimation as presented in the cited refer- 
ence. The authors concatenated 100 or 50 sets of 25 or 50 dimensional Sobol' sequence, 
respectively. They exploit the LSS method in order to obtain a complete 2500 dimen- 
sional sample of digital net. Owen 19] is more restrictive; the author suggests to use 
scrambled digital sequences for the first five or ten components and LHS for the others 
or to concatenate the principal components. We compare the results and investigate 
the effective dimensions and the contribution of the eigenvalues of the global correlation 
matrix. Table [S] reports input parameters for our test. 

We compute the eigenvalues and eigenvectors of Smw- Property (10) of the Kro- 
necker product is fundamental in this computation and considerably reduces the com- 
putational burden and time. It results that the effective dimension is 143 or 170 for 
the correlation and uncorrelation cases, respectively, are much smaller than the nom- 
inal one. Considering the first 143(170) columns, that is the first 143(170) principal 



Si{0) 


= 100 




K 


= 100 




r 


= 4% 




T 


= 1 






= 10% + ^40% 


for i = 1, . . . , 10 


Pij 


= and 40% 


for i, j = 1, . . . , 10 



Table 6: Input Parameters Used in the Second Simulation 
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Uncorrelat ion 


Standard MC 


LHS 


RQMC 


PCA 
Cholesky 


3.414(0.015) 
3.426(0.015) 


3.4546(0.0054) 
3.4323(0.0070) 


3.4438(0.0015) 
3.4518(0.0058) 


Correlation 


Standard MC 


LHS 


RQMC 


PCA 
Cholesky 


5.648(0.029) 
5.604(0.029) 


5.6655(0.0032) 
5.670(0.013) 


5.65750(0.00040) 
5.63710(0.019) 



Table 7: Prices and RMSEs both for the correlated and uncorrelated case when 100% of 
the variance is considered. 





Positive Correlation 






Zero Correlation 




Price 


RMSE 


E 


Price 


RMSE 


E 


5.262 


0.090 


5 


2.596 


0.041 


5 


5.294 


0.088 


10 


3.190 


0.047 


10 


5.433 


0.088 


15 


3.212 


0.047 


15 


5.528 


0.091 


20 


3.239 


0.047 


20 


5.484 


0.092 


25 


3.289 


0.047 


25 


5.445 


0.090 


30 


3.375 


0.048 


30 


5.653 


0.015 


147 


3.452 


0.010 


170 



Table 8: Prices and RMSE for different principal components when LHS is used. 



components, the generating matrix C takes into account 99% of the total variance. 

Table [7] shows all the results we obtained. We concatenate 50 sets of 50-dimensional 
randomized low-discrepancy sequences. 

We consider both the matrix of 2500 rows and 143(170) columns, excluding the 
effects of the remaining principal components, and the complete ANOVA in order to 
investigate the effectiveness of our assumptions and hypotheses. 

Table [S] presents the different Monte Carlo estimations with respect to the number 
of eigenvalues when LHS is used. 

Table[9]illustrates the values found by Imai and Tan [11 . Their results were obtained 
assigning the importance of each component (not anymore PCA) with their LT method. 
All the estimations found are unbiased and in agreement with those presented in the 
cited references. 

The Quasi-Monte Carlo method with LSS extension proves to be a powerful vari- 
ance reduction technique, particularly when coupled with the ANOVA-PCA decom- 
position. Moreover, the Kronecker product turns out to be a fast tool to generate 
multi-dimensional Brownian paths. Indeed, the elapsed time to realize the same path 
without using the properties of the Kronecker product is a lot higher. 

The estimation with Cholesky decomposition gives higher uncertainty than the PCA 
approach, meaning that a small amount of variance is lost. This is due to the fact that 
a relevant part of the variance is carried out by a few eigenvalues of the covariance 
matrix R. If these eigenvalues are observed, it can be noticed that only few of them are 
relevant in the PCA analysis and they are much bigger than the ones of the matrix S. 



Uncorrelat ion 


RQMC 


Correlation 


RQMC 


PCA 

Cholesky 

LT 


3.4475(0.0023) 

3.426(0.0087) 

3.4461(0.0012) 


PCA 

Cholesky 

LT 


5.65860(0.00072) 

5.603(0.022) 

5.6780(0.00047) 



Table 9: Estimated Results by Imai and Tan [13] 
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9.2 Constant Volatility: Comments 



Based on these results, we can make the foUowing conclusions: 

1. The RQMC method and the use of the Faure-Tezuka scrambling technique pro- 
vide the best estimation among all the implemented procedures for both the 
"Correlation" and "Zero Correlation" cases. The correspondent RMSEs are the 
smallest ones with a higher order of convergence with the same number of simu- 
lations. 

2. The Kronecker product is a fast and efficient tool for generating multi-dimensional 
Brownian paths with a low computational effort. 

3. As compared to to the standard Monte Carlo and LHS approaches, the use of 
scrambled low-discrepancy sequences provides more accurate results, at least for 
M X < 20, particularly with the PCA and LT-based methods. 

4. The accuracy of the estimates is strongly dependent on the choice of the Cholesky 
or the PCA approach. In particular, independent of the simulation procedure 
(MC, LHS or RQMC), when using PCA decomposition the estimates are affected 
by a smaller sampling error (smaller standard error). 



9.3 Time-dependent Volatility: Results 

The constant volatility hypothesis is the starting point for the pricing problem. A fur- 
ther improvement can be achieved by considering a time-dependent volatility function. 

It is market practice to choose step-wise time-dependent volatilities. We want to in- 
vestigate a more complex dependence to test our new approach based on the Kronecker 
product approximation. For this aim, we adopt an exponentially decaying function 
having the following expression: 

ai ^ ai{Q)exp(- —) + (Ti{+(X)) (71) 

where o'i(O) -I- (Ti(+oo) = cri(O) is the initial volatility for the i-th asset, ai{+oo) is its 
asymptotic volatility and its decay constant. 

The particular time-dependent function leads to the following solution: 



I 



(Ji{t)ak{t)pikdt = (Ti(0)(Tfe(0)Tjfe(l - ea:p( )) + 

^ Tik ^ 



+(Ti(Q)ak{+oo)n(l - exp( - — )") -|- 
+ak{0)ai{+oo)nk(l - exp( —)) + 

\ Til. / 



+a.i{+oo)ak{+oo)t 

where nk = nTk/in + Tk)- 

The simulation implemented to obtain the price of an Asian option supposing time- 
dependent volatility evolves as the constant volatility case. The main difference is the 
procedure to reduce the dimension of the problem. 

The parameters chosen for the simulation are listed in table [TOl 

The initial volatilities are equal to those used in the constant volatility case. The 
asymptotic volatility and the decay constant are the same among all the assets. These 
parameters are chosen in order to allow a comparison with respect to the constant 
volatility case. Indeed, the price of the options is sensitive to the change of volatility 
and in particular its decreasing trend should provide a lower price. 

The basket consists of 10 underlying assets, the time grid has 250 equally spaced 
points and the number of runs is = 8192 and 10 replications. Table [TT] shows the 
results coming from the simulation using the RQMC method both with the KPA and 
the PCA for dimension reduction. 
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= 100 


r 


= 4% 


T 


= lyear 




= 10% + ^40% 


(Ti( + Oo) 


= 9% for all i 


Ti 


= 1.5year 


K 


= 100 


Pij 


= and 40% for i,j = l,...,W 



Table 10: Input Parameters for the Time-depending Case 



Positive Correlation 


(KPA) 




Zero Correlation(KPA) 


Price 




5.19658 


Price 


3.20784 


RMSE 




0.00063 


RMSE 


0.00040 


E 




145 


E 


173 


Positive Correlation 


(PCA) 




Zero Correlation(PCA) 


Price 




5.19856 


Price 


3.20147 


RMSE 




0.00062 


RMSE 


0.00040 


E 




123 


E 


150 



Table 11: Estimated Results for the Time-depending Case, ANOVA = 0.99 



The KPA path-generation is efficient and fast. To have an idea of its speed, the 
elapsed times to obtain the generating matrix without exploiting the properties of 
the Kronecker product and no approximations are more than ten times higher. As 
expected, the simulation gives smaller prices with respect to the constant volatility 
situation, because a decreasing volatility function has been assigned. 

The nominal dimensions of the problem E using PCA come out to be 126 and 150 
for the correlation and uncorrelation cases. When adopting the KPA the approximated 
nominal dimensions are higher, 145 and 173, respectively. If we would consider ANOVA 
= 0.9885 for the correlation case and ANOVA — 0.98805 for the uncorrelation case 
we would get the PCA-found nominal dimension for ANOVA=0.99. We can judge 
this small difference as negligible and consequently our approximating technique to 
be efficient and leading to consistent results. As with N and M small, the Cholesky 
decomposition alone would require a small number of operations without giving any 
order of the importance for the random sources. 

Tables [12] and [T3] present the estimated prices when taking into account the full 
components. All the results are in accordance with those ones found with ANOVA= 
0.99. 

Table [TJ] illustrates the sensitivity with respect to the number of principal com- 
ponents E: As in the constant volatility case it can be seen that the estimation is 
convergent. 

Moreover, we launch a new simulation with the LHS technique with the same set 
of parameters. We list the estimated results in table [151 The estimated prices have 
higher RMSEs, confirming the fact that the RQMC approach provides a good variance 
reduction. 





KPA 


PCA 


Cholesky 


Price 
RMSE 


3.20545 
0.00040 


3.20390 
0.00041 


3.1838 
0.0091 



Table 12: Uncorrelation Case. Estimated Prices and Standard Errors. ANOVA= 1. 



30 





KPA 


PCA 


Cholesky 


Price 
RMSE 


5.20060 
0.00050 


5.20210 
0.00058 


5.1946 
0.0093 



Table 13: Correlation Case. Estimated Prices and Standard Errors. ANOVA= 1. 





Positive Correlation 






Zero Correlation 




Price 


RMSE 


E 


Price 


RMSE 


E 


5.7805 


0.0079 


5 


2.6368 


0.0038 


5 


4.9904 


0.0081 


10 


2.9681 


0.0042 


10 


5.0226 


0.0081 


15 


3.1172 


0.0043 


15 


5.1103 


0.0081 


20 


3.0979 


0.0043 


20 


5.1826 


0.0083 


25 


3.1051 


0.0043 


25 


5.1937 


0.0082 


30 


3.1514 


0.0043 


30 



Table 14: Prices and RMSEs for different principal components. Case: RQMC 





Positive Correlation 






Zero Correlation 




Price 


RMSE 


E 


Price 


RMSE 


E 


4.874 


0.016 


5 


3.121 


0.089 


5 


5.093 


0.016 


10 


3.118 


0.085 


10 


5.097 


0.016 


15 


3.122 


0.086 


15 


5.131 


0.016 


20 


3.072 


0.086 


20 


5.145 


0.016 


25 


3.163 


0.088 


25 


5.201 


0.016 


30 


3.110 


0.089 


30 



Table 15: Prices and RMSEs for different principal components. Case: LHS 
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9.4 Time-Dependent Volatility: Comments 



According to the results we have found in the time-dependent case, it is possible to 

draw the following conclusions: 

1. RQMC with LSS is a general approach that does not depend on the chosen price 
dynamic. 

2. The KPA we propose, provides unbiased estimations with a reduction of the 

computational cost. In the framework we investigate, KPA returns a higher 
nominal dimension, as expected, but only relatively to a negligible amount of 
variance. 

3. KPA is a lot faster than the straightforward PCA because it exploits the prop- 
erties of the Kronecker product and the boomerang shaped matrices. The ad hoc 
Cholesky decomposition algorithm we develop is fundamental for the KPA. We 
do not report computational times because we expect that further improvements 
can be done. 

4. KPA and PCA can be considered both valid as path-generation methods to sup- 
port the ANOVA and the identifications of effective dimensions. 
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